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Abstract 

A mean field analysis of finite density QCD is presented including the effects 
of additional chiral invariant four-fermion interactions. A lattice regulariza- 
tion is used with Nj = 4 flavors of staggered fermions. The use of the four- 
fermion coupling as an improved extrapolation parameter over the bare quark 
mass in Monte Carlo simulations is discussed. Particular attention is given 
to the structure of the phase diagram and the order of the chiral phase tran- 
sition. At zero gauge coupling, the model reduces to a Nambu-Jona-Lasinio 
model. In this limit the chiral phase transition is found to be second-order 
near the zero-density critical point and otherwise first-order. In the strong 
gauge coupling limit a first-order chiral phase transition is found. In this 
limit the additional four-fermion interactions do not qualitatively change the 
physics. The results agree with previous studies of QCD as the four-fermion 
coupling vanishes. 
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I. INTRODUCTION 



Monte Carlo simulations of lattice gauge theories with dynamical fermions become very 
computationally intensive at small fermion masses due to the zero mass singularity of the 
Dirac operator. In order to study theories with light fermions it becomes necessary to carry 
out simulations with relatively large bare masses, and then extrapolate back to lower values. 
In studying models with chiral symmetry or dynamical chiral symmetry breaking it would 
be desirable to simulate the models directly in the chiral limit, without having to introduce 
any bare mass terms which explicitly break the symmetry. 

A new approach to this problem has been suggested f^-^ which does not explicitly break 
chiral symmetry. The idea is to add additional terms to the lattice action, corresponding 
to chiral invariant four-fermion interactions. In the continuum limit, as the lattice spacing 
reduces to zero, these additional terms are irrelevant (in the renormalization group sense). 
However, at finite lattice spacings they favor the generation of a dynamical fermion mass. 
This dynamical mass removes the singularity in the Dirac operator allowing for simulations 
with exactly zero bare fermion masses 

One particular area where this may soon lead to significant progress is in the study of 
QCD at finite density. It is generally expected that at some density greater than that of 
normal nuclear matter there is a transition to a quark-gluon plasma and chiral symmetry is 
restored. The nature and order of this transition may be an important factor in cosmological 
models, neutron stars, and future heavy ion collision experiments. Mean field calculations 
at both infinite and strong gauge coupling indicate that the transition is first-order 
at low temperatures. Unfortunately, Monte Carlo simulations of QCD at non-zero density 
have not been able to give a physically accurate description of the transition Q,^]. 

Quenched simulations of QCD using the grand canonical formulation show the tran- 
sition at a value of the chemical potential equal to one half the pion mass, vanishing as 
the bare quark mass mo 0. These results are now understood to correspond to the 
Nf^O limit of a theory of Nf quarks with regular action and Nf quarks with the conjugate 
action, and are therefore unrelated to QCD. 

Unquenched Monte Carlo calculations have to deal with a complex fermionic determi- 
nant due to the introduction of the chemical potential. Hence, conventional Monte Carlo 
algorithms can not be applied. Some results have been obtained at infinite gauge coupling 
by representing the partition function in terms of monomers, dimers, and baryonic loops 



||1 1|| , and also on small lattices, by using a spectral density method |T2[. The results are 
consistent with mean field calculations. However, it is not clear how to extend these methods 
to intermediate or weak gauge coupling. 

The evaluation of the complex determinant can also be avoided by the fugacity expansion 
method of Refs. [p!3l , |T4| . However, the proper physics of chiral symmetry breaking, with a 
massless pion, has not yet been achieved using finite quark masses with this method [0. 
Simulations of zero mass QCD with additional chiral four-fermion interactions, called xQCD 
1^, may be able to overcome this problem. In support of this approach, this paper presents 
a mean field analysis of xQCD at non-zero density. The calculations will be restricted to the 
two limiting cases of zero and infinite gauge coupling. No restrictions are placed the strength 
of the four fermion interactions. Although in applications approximating QCD they should 
be chosen relatively weak. 
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II. THE xQCD MODEL 



The xQCD action, using the staggered fermion |T^ formulation, is 

S = Sf + Sg, 



X 

Nf/4 

a=l '-x,y X <x,x> 

= -/5E[1 - ^Re(Trnf/f/t/t/)]. (1) 

Xi and Xi are complex Grassmannian fields defined on each site of the lattice. Each Xi 
represents four degenerate fiavors of fermions, so the sum over the fiavor index a runs from 
1 to Nf/A. We are interested in studying QCD, however, the calculations can be done for 
any SU(A^) gauge group. Considering this more general case, the color indices i,j run from 
1 to N. 

The fermion hopping matrix Ai^^{x,y) is given by 



M'^ix,y) = ^j2[U{x)Ui^{x)6y,,^, - f-\x)Ui^{x - 9)6y^,.o] +mo^,,., (2) 

with 

- I T^C^ ' ^ ? 9 . • (3) 

[ Vu{x) u = 1,2,3 ^ ^ 

Here fi is the chemical potential [p!6|,|6[] measured in units of the lattice spacing. The T]y{x) 
are the conventional Kawamoto-Smit phases (^—i"j^o+---+x^-i ^ jg ^^iq bare (current) 
mass of the quarks which can be set to zero, but is included here for completeness. 

The auxiliary fields a and vr have been introduced to separate the four-fermion inter- 
actions [|T^,0. They are defined on the dual lattice. The symbol <x,x> represents the 
16 dual sites x adjacent to the direct lattice site x. e{x) represents the alternating phase 
(■_]^^a;o+xi+x2+a;3_ ^f^gj^ integrating over the auxiliary fields, the lattice model approximates 
the continuum theory 

C = i,{]p + /i7o + mo)^ - 2]^[(^^)' - i^l^m - ^TiiF.^F'^'^), (4) 

with G = l/'~f and 2N/g'^ = (3. (The spinor, color and fiavor indices have been suppressed.) 

In the chiral limit (mo = 0), the continuum theory (^ is invariant under the global U(l) 
chiral symmetry 

In the language of the lattice model this translates to 
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cos(2^) sm{2e) \fa\ , . 

-sm{29) cos(2e) )\7r ) ' ^ ' 

The U(l) chiral symmetry is spontaneously broken in the low density hadron phase through 
the dynamics of the color gauge fields. The chiral condensate 

(XV) = ^7(cr), (7) 

is nonzero in the broken phase, and serves as the order parameter for the chiral phase 
transition. 



and 




III. ZERO GAUGE COUPLING NJL LIMIT 

As the gauge coupling vanishes xQCD reduces to a lattice Nambu-Jona-Lasinio (NJL) 
model . In Ref. ||2l| the three dimensional version of this model was studied at non-zero 



chemical potential in a 1/Nf expansion. Presented here are results for four dimensions. Of 
particular interest is the structure of the finite density phase diagram and the order of the 
chiral phase transition. 

When f3 —y oo the SU(A^) gauge variables are frozen to Ul^{x) = 6^^ at each link of the 
lattice. Gaussian integration over the fermion fields leaves the following, bosonic, partition 
function 

Z = l DaDvr exp|-^7E(^'(5) + ^'(5)) + ^Tr ln[5^i]}, (8) 



where 

3 

2 



1 ^ 

Sp^i^^y) = [fAx)^y,x+i> - f^\x)6y,^.o] +mo+j^ J2 {(^{x) + ie{x)7T{x)). (9) 

u=0 <x,x> 

The mean field, saddle point solution {a, n} (equivalent to the first order of a large Nf 
expansion) is given by the "gap equation" 

(^)7=('*V"°)™'|5-1- (1") 
where (in the infinite volume limit) 

d^A; 1 



Trf^i?] = 16 - — -r- = ^ 

J-f (27r)4 i(l - cos 2A;o cosh 2/i - i sin 2A;osinh2/i) + 22 sin^ + 

v=l,2,3 

,2 I _ \2 I -2 



2 ' 

m 



= {a + moY + r. (11) 

Here m is the effective (dynamical) mass of the quarks. If the bare quark mass mo is non- 
zero, chiral symmetry is broken explicitly. The unique solution is then identified with a>0 
and 7r = 0. In the chiral limit the solution is invariant under the transformation (|^), and so 
we are free to choose 7f = 0. This is the case of interest, which we will now discuss for N = 3. 
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The gap equation is easily evaluated numerically on a finite size lattice by replacing the 

A (48)^ Euclidean lattice was found to 



momentum integrations with discrete sums ^1 



be sufficiently large to eliminate finite-size effects. Figure |I] shows the resulting behavior 
of a at zero chemical potential. The other data points show (zero mass) measurements of 
+ from Monte Carlo simulations of xQCD in the P^oo limit ||2^. The mean field 
phase transition is at 7c = 1.86. This is also in agreement with Monte Carlo simulations of 
the pure NJL model p3, where the transition was found at 7c = 1.81 {'Jc/N jSc in their 
notation). 

At non-zero chemical potential the chiral phase transition is given by setting a = 0. The 
gap equation then produces the critical line shown in the phase diagram of Fig. |^. The 
upper left section of the phase diagram corresponds to the spontaneously broken symmetry 
phase, with a > 0. The remainder of the diagram represents the chiral symmetric phase, 
with a = 0. 

As an alternative to solving ([T0|) the minimum energy solution can be found by a direct 
evaluation of the action. This also allows for a clearer determination of the order of the 
chiral phase transition. Substituting the uniform ansatz {a, 0} into the action and Fourier 
transforming one finds 



S= 17^2-3-16 



In 



1 — cos 2ko cosh 2/i — i sin 2kQ sinh 2/i) 



sin^ k^ 



a 



=1,2.3 



(12) 



Figure § contains plots of S* as a function of a at 7 = 1. Each line corresponds to a different 
value of the chemical potential. The plots have been rescaled so that 5* is zero at a = 0. 
The figure shows a clear first-order (discontinuous) transition as the chemical potential is 
increased from zero. However, at larger values of 7 the transition becomes continuous, as 
shown in Fig. |^ for 7 = 1.5. The chiral transition remains continuous as the critical line is 
followed towards /z = 0. This appears consistent with the second-order transition found in 
the 1/A^/ expansion at zero-density ||2^.^ What is perhaps unanticipated is that the second- 
order transition actually extends into the interior of the phase diagram for small values of 
the chemical potential. 

The behavior of a as a function of the chemical potential is plotted in Fig. ^, for the 
fixed values of 7 = 1 and 1.5. Other quantities of interest include the quark number density 
and quark energy density. 



d\nZ 



-NTr 



OfJ, 



^ Monte Carlo calculations |23| indicate that the theory is in fact trivial, exhibiting logarithmic 
scaling corrections to mean field exponents. 
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= 3-16 



f d^k |(cos2fcosinh2yU + isin2fcocosh2yu) 

f (27r)4 1(1 _ cos 2A;o cosh 2/i - i sin 2/co sinh 2/i) + ^ sin^ A;,, + ct^ ' 



i^^doij) = NTt[{Sp')'Sf] 

rf d'^k i(l — cos 2A;o cosh 2/i — i sin 2A;o sinh 2/i) 

= 3-16 / - — -r- = T, t:- (13) 

J-f (27r)'' i(l - cos 2A;o cosh 2/i - i sin 2A;o sinh 2/i) + 22 sin^ + ^ 

i/=l,2,3 

These are plotted in Fig. ^ and Fig. 0, respectively, for the values of 7 = 1 and 1.5. As 
expected, the quark density is essentially zero in the broken phase. As the chemical potential 
is increased towards the critical value separating the two phases, roughly equal to m = a, the 
order parameter begins to decrease and the density rises. At the chiral transition the quark 
mass abruptly vanishes. Hence, the quark density abruptly increases. In the symmetric 
phase, the number density increases ~ /i^ at low densities, which is consistent with a non- 
interacting quark gas. As the lattice becomes filled the density eventually saturates at the 
maximum value allowed by the Pauli exclusion principle, of A^ = 3. 



IV. THE STRONG GAUGE COUPLING LIMIT 

In the previous section we studied the case where the gauge coupling vanished. However, 
in using xQCD to approximate QCD one is really interested in the case where the four- 
fermion interactions are very weak compared to the gauge interactions; naively when G<^ g. 
In this section we will study the strong gauge coupling limit of g—>-oo. Similar infinite gauge 
coupling calculations have been carried out in [^, as well as in [|],|| at finite density. Some 
of the same techniques will be used here with suitable modifications. 

Begin by defining the following "meson" and "baryon" fields: 

M^'ix) ^ xr(x)x^(x), 

B'^'-^x) ^ l^e,,...,x-{^)x%x) ...xlix). (14) 

The meson fields commute for any value of N, while the baryon fields will anticommute 
(commute) for odd (even). 

At infinite gauge coupling (/? = 0) the gauge term in the action Sg is zero. Therefore, the 
gauge variables on each link are completely independent. Integrating over the SU(A^) gauge 
group at each link yields^ 

DxDxDaDTT 



^ Equation (|^) is valid only for A^>3. Similar expressions for the U(l) and SU(2) gauge groups 
can be found in Ref. p4 
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r N ^^^^ r 

X expj- ^7 5^(^2 (£) + 7r2(x) )- Ex"(a;)x-(x) + ^ ^ + i£(x)7r(£)) 

^ i a=l X L <x,x> 

+ 32iV(Ar2-l) (Tr/[M(x)M(x + z>)]^ + iVTr^[M(x)M(x+z>)M(x)M(x + z>)]) + . . . 



'_i\Ar(Ar+l)/2 



• (15) 



Here Tr/ is the trace over flavor indices only. In writing (|T^), the full expression for the 
meson contributions has been truncated. (In general the series will terminate after terms 
of order 0(M^^^/^) due to the Grassmannian nature of the Xi fields.) For the remainder 
of the calculation we only need to retain the first meson-meson term. As shown in 



this is consistent with a systematic expansion in l/rf, where d is the number of space-time 
dimensions. 

Exponentiating the product at each lattice site yields the following action: 



'5 = ^7E(^'(5)+^'(5^)) + ETr/[M(x)] + E {<y{i) + H^)<x)) 

X X ^ <x,x> ' 

- Ef|Tr;[^(a^)^(^'2/)^(l/)] + Tr/[5(x)FB(x,2/)i?(2/)]), (16) 



x,y 



with 



1 3 

y) = jrT E( + ^y,x-u), (17) 

/ i\JV(JV+l)/2 3 

VB{x,y) ^ ^—h^ E(/.^(^) W + (18) 

Note that all the dependence on the chemical potential is now contained in VB{x^y). 

We want to construct an effective action by integrating out the staggered fermion fields, 
as was done in the pure NJL case. To facilitate this we first linearize the action in the meson 
and baryon terms using the following identities: 



x,y 



x,y 



exp^EiTr/[M(x)F(x,2/)M(2/)] 

det[V] /Dp e^J-Y,\TTf[p{x)V-\x,y)p{y)]-Y.p'^\x)M'^\x)\, (19) 
exp^^Tr^[5(x)FB(x,i/)i?(y)] 

dei[VB\ him e^J-Y,i:if[h{x)V^\x,y)h{y)] 

^,y 

X expj- ^(6"''-^(x)5"^-"(x) + 5"^-"(x)6"*-"(a;))|. (20) 
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In writing the second identity we have taken to be odd, as it is for QCD. In either case, 
the final expression for the effective action (|2^) will be correct for all iV > 3.0 

Let us now restrict the remainder of the calculation to Nf = 4 flavors which is techni- 
cally much easier to handle. (Unlike the high temperature transition (at low density), the 
finite density transition should not have a large flavor dependence at low temperatures |0].) 
Integrating over xt and Xi we then flnd 

Z = jDaDnDpmm exp|-^7 ^(^^(5) + 7r2(£))| det[y] det[yB] 



X det 



p(x) +mo + (^(^) + i£{x)n{x)) 



<x,x> 



N 



N{N-1) - ^ 

-1) 2 b{x)b{x) 



X 



x,y 



(21) 



Since det[y] does not depend upon any of the flelds or the chemical potential, it can be 
dropped without affecting the calculation of any physical quantities. The remaining func- 
tional determinants can be exponentiated up into the action, producing an effective action 



x,y 



Seff = hY.i^'{x) + 7r\x)) + Y.[lpix)V~\x,y)p{y)+b{x)V^\x,y)b{y) 



Trln[S^(x) + {-l)^^^bb] - TrlnfFe] 



(22) 



with 



E(x) = p(x) -|- mo + 1^ X! ('^(^) + i£{x)7!'{x)). 



(23) 



<x,x> 



Using the fact that ln[l + ^] = ^ for the Grassmannian flelds b and b, ( p2D can be further 
simplifled to 



Seff = llT.i^\x) + n\x))+J2np{x)V-''ix,y)piy) + bix)S^\x,y)biy) 



Trln[S^(x)] -Trln[yB], 



(24) 



with 



S;\x,y) = V^\x,y) + ^^. 



(25) 



Integrating over b and b, which are now just Gaussian integrals, we arrive at a flnal expression 
for the effective action at /3 = 



^ For even det[yB] is replaced by det[yg^] in (pO[). Integration over 6,6 in ( [2l| ) then gives 
J LhLbbbe~^^B * = det[yg]. Combining all the determinants eventually yields equation (^ 
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Seff = hE(^'(^) + ^'(^)) + E lp{x)V-\x,y)piy) - Trln[S^(x) + VB{x,y)]. (26) 



x,y 



Now consider a uniform (mean field) saddle point approximation to the action. Denote 
the saddle point solution which minimizes the effective action by {p, a, vf}. Substituting this 
into (p6|) and transforming to momentum space yields 



In 



cos 2ko cosh 2iV/i — i sin 2ko sinh 2NiJ,) 

(27) 



+ ^ sin^ + m| 

v=l,2,3 



with 



m 



iN-lu- 



B 



[{p + a 



(28) 



Seff can be evaluated by an explicit summation over discrete momenta on a finite sized 
lattice, as carried out in section |T1|. In the chiral limit (mo = 0) we find two distinct phases: 
For pKpcipi) there is a spontaneously broken phase withQ 



P 



N ' 



N _ 



TT 



0; 



while for fi> fidl) there is a chiral symmetric phase, characterized by 

p = a = 71 = 0. 



(29) 



(30) 



(The same solutions are produced by evaluating the "gap equation" , obtained by minimizing 
(P7|) with respect to {p, a, 7f }.) Notice that within each phase all of the order parameters are 
independent of the value of the chemical potential. Therefore, the chiral phase transition is 
discontinuous (first-order) for all values of 7. The critical line can be determined by setting 
Seffiip, ^> Pc) = SeffiiO, 0, 0}; iXc)- This yields 



/i,(7) = iln[4(2 + iV/7)]-l 



2' 



(31) 



in agreement with previous mean field QCD calculations P,^] as 7 — > 00. 

Let us now discuss the specific gauge group SU(3). The phase diagram is shown in Fig. ||. 
As expected, the values of pc increase monotonically as G = l/7 increases. Of course, since 
g^G, chiral symmetry breaking is still due to the strong gauge (color) forces as in ordinary 
QCD. The additional four-fermion interactions merely enhance the effect. (This is clearly 
demonstrated by comparing Fig. ^ with the phase diagram produced in the absence of the 
gauge fields (Fig. 0).) The mean field chiral condensate is given by 



^ Since the values of {p, a, vf} are independent of ft in each phase the analytic expressions ( p9D 
can be obtained with fi = 0. In which case, the baryon terms can be safely ignored and VB{x,y) 
can be set equal to zero. 



9 



(32) 



Values of p and a in the broken phase are plotted as a function of G= I/7 in Fig. |^. Included 
for comparison are values of a from exploratory Monte Carlo simulations of xQCD at 
/3 = 0.5 and zero chemical potential. 

Plots of the effective action for 7 — >■ cx) are shown in Fig. [10| at selected values of the 
chemical potential /i. It should be obvious from the figure that the chiral phase transition 
is indeed discontinuous, and that p does not depend on the value of the chemical potential 
within the broken phase. The value of /ic = 0.54 agrees with the mean field result of for 
ordinary QCD. Figure |Tl] contains similar plots of the effective action for 7 = 1. Notice that 
the general behavior of the action and the order of the transition remain unchanged. (This 
was explicitly checked for values down to 7 = .01.) 

Let us now consider the mass of the baryons in the broken phase. By adding baryon 
source terms to the original action, it is straightforward to show that the staggered baryon 



propagator {B{x)B{y)) has the same pole structure as Sb{x,y), defined in (|25|). Defining 
the baryon mass Mb{p) to be the location of the pole of {B{x)B{y)) in the broken phase 
gives 



Mb{p) = sinh ^(m^) — Np = sinh 



4(2 + f 



3p. 



(33) 



In the symmetric phase the dynamical quark mass vanishes along with the chiral order 
parameter. The pole in {B{x)B{y)) disappears, and so, as expected we no longer have 
stable nuclear matter. 

Mb{p) is plotted in Fig. |T2| for various values of 7. Notice that the chiral symmetry 
restoration occurs consistently at a value of Mb{pc) = 3/2. This value can be traced back 
to the first two terms of Seff in (pTf), which represents the mean field energy density of 
the chiral condensate. At weaker gauge coupling, where the binding energy is small, one 
would expect the phase transition to be closer to the baryon mass threshold Mb{pc) = 
0]. In ordinary QCD, one finds that this energy difference does indeed decrease when l/g"^ 
corrections are added to the calculations [Q. 

The quark number density and quark energy density are given by the following equations: 



-Tr 



3-16 



+ Vb 



I (cos 2ko sinh 2Np + i sin 2 fco cosh 2Np) 



{ij^tP) = Tr 



'-I (27r)4 1(1 _ cos 2ko cosh 2Np - i sin 2ko sinh 2Np) + ^ sin^ k^ + 



m 



.2 ' 

B 



16 



f 

-f (27r)'^ 1(1 - cos 2fco cosh 2A^/i - i sin 2A;o sinh 2A^/i) + ^ sin^ fc,, + 

i^=l,2,3 



1(1 — COS 2/co cosh 2Np — i sin 2kQ sinh 2Np) 



In the broken phase the baryon mass gap remains relatively large all the way to the phase 
boundary. Therefore, the baryon density is strictly zero. In the symmetric phase the binding 
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energy of the quark matter is very large at strong coupling. Hence, the density jumps directly 
to saturation as the chemical potential is increased through the transition. This is certainly 
an artifact of working at infinite gauge coupling. Similar behavior has been found in other 
infinite coupling calculations 0. 
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FIG. 1. The solid line shows a vs. G= I/7 as given by the /3^oo gap equation at zero chemical 
potential. The data points are Monte Carlo measurements of Vc^ + vr^ at zero mass from Ref. [22|. 
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FIG. 2. The mean field phase diagram of xQCD aX (3 ^ 00. Plotted here is the critical line 
separating the dynamically broken phase, at large G, from the chiral symmetric phase, at small G. 
G = l/7. 
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FIG. 3. The (3 ^ oo effective action of xQCD plotted as a function of a, for 7 = 1.5. The 
different hnes are for /x = 0, /x = .60, /x = .69, and /x = .72, from bottom to top. 
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FIG. 4. The (3 ^00 effective action of xQCD plotted as a function of a, for 7= 1. The three 
different lines are for /x = 0, /x = .90, /x = .97, and /x=1.0 from bottom to top. 
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FIG. 5. cr as a function of the chemical potential at /3 — >^ oo for 7 = 1.5 (solid line) and 7=1 
(dashed line). 
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FIG. 6. Plots of the quark density at /3 — 00 for 7 = 1.5 (solid line) and 7 = 1 (dashed line), as 
a function of the chemical potential. 
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FIG. 7. Plots of the quark energy density at /3 — oo for 7 = 1.5 (solid line) and 7 = 1 (dashed 
line), as a function of the chemical potential. 
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FIG. 8. The mean field phase diagram of xQCD at /? = 0. The solid line is the critical line 
separating the dynamically broken phase at small from the chiral symmetric phase at large fj,. 
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FIG. 9. Mean field values of p (solid line) and a (dashed line) as a function of G = I/7 in the 
broken phase, at /3 = 0. The diamonds are Monte Carlo values of a at /3 = 0.5 from Ref. [22[|. 
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FIG. 10. The /3 = effective action of xQCD plotted as a function of for 7— >oo. The three 
different lines are for /i = 0, /i=.54 and /i = .7, from bottom to top. 
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FIG. 11. The (3 = effective action of xQCD plotted as a function of p, at 7 = 1. The three 
different hnes are for /x = 0, /i = 1.0 and /x = 1.4, from bottom to top. 
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FIG. 12. Plots of the (3 = baryon mass as a function of the chemical potential in the broken 
phase. The three lines are for 7— >-oo, 7 = 10, and 7 = 1, from left to right. 
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